home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / mahe / ahecalc_short.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-17  |  10.6 KB  |  404 lines
  1. /* 5/30/1988 Modified by Shie-Jue Lee:                */
  2. /*    1. data type for equal_const                */
  3. /*    2. Iteration controls for histogram calculation        */
  4.  
  5. #include <stdio.h>
  6. #include <signal.h>
  7. #include <math.h>
  8. #include "ahe.h"
  9.  
  10. #define Malloc(x) malloc((unsigned)(x))
  11.  
  12. #define iptr1(Y,X)    (im1+((Y) * numline)+(X))
  13. #define iptr2(Y,X)    (im2+((Y) * numline)+(X))
  14.  
  15. #define mapval(J,K)    (*(maps + ((J) * regline) + (K)))
  16.  
  17. int         status;        /* progress of calc. indicators     */
  18.  
  19. ahecalc_short(im1, im2, dimv, minmax, nreg, clipfrac, argv)
  20.      /* CHANGED */
  21.      unsigned short *im1;            /* input image     */
  22.      unsigned char  *im2;            /* output image */
  23.      int *dimv,            /* array of image bounds     */
  24.      minmax[2],        /* minimum and maximum inputs     */
  25.      nreg[DIMNO];        /* number of regions in x, y, z  */
  26.      float clipfrac;        /* fraction of average bucket cnt */
  27.      char **argv;        /* command line for status report via 'ps' */
  28. {
  29.  
  30.  
  31. /* CHANGED */
  32.     int tmpint ;
  33.     int         min,        /* minimum grey value             */
  34.                 max,        /* maximum grey value             */
  35.                 maxlevels,    /* number of grey levels in image     */
  36.                 numregions,    /* number of regions to be used         */
  37.                *del[DIMNO],    /* region sizes in each dimension     */
  38.                *curdel,        /* current region under consideration     */
  39.                 regsize,    /* compute size of current region     */
  40.                 excess,        /* excess pixels in each dimension     */
  41.                 newdel[DIMNO],    /* size of current regions         */
  42.                 start[DIMNO],    /* current starting position in image     */
  43.                 regline,    /* number of regions in line of image     */
  44.                 cliplimit,    /* computed limit of bucket count     */
  45.                 limit,        /* cliplimit - ptemp             */
  46.                 numline,    /* number of pixels per line of image     */
  47.                *histo,        /* histogram array             */
  48.                 ptemp;        /* floor of p                  */
  49.     float      *cumhisto,    /* cummulative histogram array           */
  50.                 p;        /* base for cumulative histogram     */
  51.  
  52.     /* CHANGED */
  53.     register unsigned short *i1;    /* fast pointers to subsegments of      */
  54.     register unsigned char *i2;        /* image arrays */
  55.  
  56.     int         ymin,        /* bounds for mapping of image         */
  57.                 ymax,
  58.                 xmin,
  59.                 xmax;
  60.  
  61.     int         jval,        /* map indicies                 */
  62.                 kval;
  63.     double      recy,        /* constants to step thru image in mapping */
  64.                 recx,
  65.                 dy,
  66.                 dx,
  67.                 cdy,
  68.                 lefty,
  69.                 righty;
  70.  
  71.     float      *m1,        /* four maps for bilinear interpolation     */
  72.                *m2,
  73.                *m3,
  74.                *m4;
  75.  
  76.     float     **maps;        /* pointers to intensity maps         */
  77.     float       equal_const;    /* constant to produce histogram equal. */
  78.  
  79.     int         i,
  80.                 j,
  81.                 k,
  82.                 loop;        /* utility infielders         */
  83.     int         x,
  84.                 y;        /* ditto             */
  85.  
  86.     register int *h;        /* fast pointer to histogram     */
  87.     register float *c;        /* fast pointer to cum. histo.   */
  88.     register float *m;        /* fast pointer to current     */
  89.                     /* intensity mapping         */
  90. #ifdef godot
  91.     extern int  report ();    /* status reporting routine     */
  92.     signal (SIGUSR1, report);    /* attach report routine to sig. */
  93. #endif
  94.  
  95. /* Preliminary setup    */
  96.     status = 0;
  97.     min = minmax[0];
  98.     max = minmax[1];
  99.  
  100. #ifdef DEBUG
  101.     fprintf (stderr,"ahe: min = %d\n", min);    /***/
  102.     fprintf (stderr,"ahe: max = %d\n", max);    /***/
  103. #endif
  104.  
  105.     maxlevels = max - min + 1;    /* runs from 0 to max + 1 */
  106.  
  107. #ifdef DEBUG
  108.      /*###*/ fprintf (stderr, "maxlevels = %d\n", maxlevels);
  109. #endif
  110.  
  111. /* calculate constants for matrix indexing    */
  112.     numline = dimv[2];        /* pixels in a line     */
  113.     regline = nreg[2];        /* regions in a line     */
  114.  
  115. /* Find the total number of regions to be considered:    */
  116.     numregions = 1;
  117.     for (i = 0; i < DIMNO; i++)
  118.     numregions *= nreg[i];
  119.  
  120. #ifdef DEBUG
  121.      /*###*/ fprintf (stderr, "maxlevels = %d\n", maxlevels) ;
  122.      /*###*/ fprintf (stderr, "numline = %d\n", numline);
  123.      /*###*/ fprintf (stderr, "regline = %d\n", regline);
  124.      /*###*/ fprintf (stderr, "numregions = %d\n", numregions);
  125. #endif
  126.  
  127. /* Allocate the histogram map.                    */
  128.     histo = (int *) Malloc(sizeof (int) * maxlevels);
  129.  
  130. /* Allocate the cummulative histogram map.            */
  131.     cumhisto = (float *) Malloc (sizeof (float) * maxlevels);
  132.  
  133. /* Allocate the intensity mapping tables:            */
  134.     maps = (float **) Malloc (sizeof (float *) * numregions);
  135.  
  136.     for (i = 0; i < numregions; i++)
  137.     *(maps + i) = (float *) Malloc (sizeof (float) * maxlevels);
  138.  
  139. /* 
  140.  * Compute the mappings for each region.  We step through
  141.  * the regions, calulating the histogram and the histogram
  142.  * equalizing mapping for each region.  The same histogram
  143.  * array is used over for each region, but each region has
  144.  * its own map. 
  145.  */
  146.  
  147. /* 
  148.  * Calculate region sizes.  If the image is not evenly divisible,
  149.  * the excess pixels are distributed among the first dimv mod n
  150.  * regions
  151.  */
  152.  
  153. /* Allocate the arrays to hold the size of each region:     */
  154.     for (i = 0; i < DIMNO; i++)
  155.     {
  156.     del[i] = (int *) Malloc (sizeof (int *) * nreg[i]);
  157.     curdel = del[i];
  158.     regsize = dimv[i] / nreg[i];
  159.     excess = dimv[i] % nreg[i];
  160.     for (j = 0; j < nreg[i]; j++)
  161.         *(curdel + j) = (j < excess) ? regsize + 1 : regsize;
  162.     }
  163.  
  164. /*
  165.  *    Determine cliplimit from range of intensities in image and
  166.  *    the region size.  Image size/(number intensity levels * number
  167.  *    of regions) gives approximation of number of pixels at each
  168.  *    intensity in a region.  Clipfrac is used to scale the result.
  169.  */
  170. #ifdef DEBUG
  171.     fprintf (stderr, "ahe: clipfrac %f\n", clipfrac);
  172. #endif
  173.     if (clipfrac > 0.)
  174.     {
  175.     cliplimit = clipfrac * (dimv[0] * dimv[1] * dimv[2]) /
  176.         (maxlevels * numregions) +.5;
  177. #ifdef DEBUG
  178.     fprintf (stderr, "ahe: cliplimit %d\n", cliplimit);
  179. #endif
  180.     }
  181.     else
  182.     /* 0 is flag for no clipping */
  183.     cliplimit = 0;
  184.  
  185. #ifdef DEBUG
  186.     fprintf (stderr, "ahecalc: clipfrac %f cliplimit %d\n", clipfrac, cliplimit);
  187. #endif
  188.  
  189. /* Set up to calculate the histogram and mappings.  */
  190.  
  191.     status = 0;
  192.     start[1] = start[2] = 0;
  193.     for (j = 0; j < nreg[1]; j++)    /* loop thru y regions     */
  194.     {
  195.         for (k = 0; k < nreg[2]; k++)    /* loop thru x regions     */
  196.         {
  197. #ifdef REPORT
  198.         printf ("status: %d/192", status);
  199.         fflush (stdout);
  200. #endif
  201.         sprintf (argv[0], "status:%d/192", status++);
  202.         newdel[1] = *(del[1] + j);    /* current region in y     */
  203.         newdel[2] = *(del[2] + k);    /* current region in x     */
  204.         equal_const = (double) (maxlevels - 1) / (newdel[1] * newdel[2]);
  205.  
  206. #ifdef DEBUG1
  207.          /*###*/ fprintf (stderr, " equal_const %f\n", equal_const);
  208.          /*###*/ fprintf (stderr, " newdel[1] %d\n", newdel[1]);
  209.          /*###*/ fprintf (stderr, " newdel[2] %d\n", newdel[2]);
  210. #endif
  211.  
  212.         /* zero histogram array                 */
  213.         for (h = histo; h < histo + maxlevels; h++)
  214.             *h = 0;
  215.         h = histo;
  216.  
  217.         /* calculate histogram for this region         */
  218.             for (y = start[1]; y < start[1] + newdel[1]; y++)
  219.             for (x = start[2]; x < start[2] + newdel[2]; x++)
  220.                 ++h[*(iptr1 (y, x))];
  221.  
  222.         /* clip histogram if user requested a limit         */
  223.         if (cliplimit > 0)
  224.         {
  225.             aheclip (histo, maxlevels, cliplimit, newdel[2], newdel[1],
  226.                  &ptemp, &p);
  227.  
  228.             /* calculate cumulative histogram             */
  229.             limit = cliplimit - ptemp;
  230.             if (*histo == limit)
  231.             *cumhisto = (float) cliplimit;
  232.             else
  233.             *cumhisto = (float) *histo + p;
  234.  
  235.             for (h = histo + 1, c = cumhisto + 1; h < histo + maxlevels;
  236.              h++, c++)
  237.             {
  238.             if (*h == limit)
  239.             {
  240.                 *c = *(c - 1) + (float) cliplimit;
  241.             }
  242.             else
  243.             {
  244.                 *c = *(c - 1) + (float) *h + p;
  245.             }
  246.             }
  247.         }
  248.         else
  249.         {
  250.             /* no clipping requested */
  251.             *cumhisto = (float) *h;
  252.             for (h = histo + 1, c = cumhisto + 1; h < histo + maxlevels;
  253.              h++, c++)
  254.             {
  255.             *c = *(c - 1) + (float) *h;
  256.  
  257.             }
  258.         }
  259.  
  260. #ifdef DEBUG1
  261.          /*###*/ fprintf (stderr, "j %d, k %d\n", j, k);
  262. #endif
  263.  
  264.         /* calculate intensity mapping number         */
  265.         m = mapval (j, k);
  266.  
  267.         *m = 0;
  268.         for (c = cumhisto; c < cumhisto + maxlevels; ++c, ++m)
  269.             *m = equal_const * (*c);
  270.  
  271.         /* increment starting point                 */
  272.         start[2] += newdel[2];
  273.         }
  274.         start[2] = 0;
  275.         start[1] += newdel[1];
  276.     }
  277.  
  278. /* Apply mappings to the image */
  279.     start[0] = 0;
  280.     loop = 0;
  281.  
  282.         start[1] = 0;
  283.         ymin = 0;
  284.  
  285.         for (j = -1; j < nreg[1]; j++)    /* y regions */
  286.         {
  287.         if (j == -1)
  288.         {
  289.             jval = 0;
  290.             dy = 0;
  291.             recy = 0;
  292.             ymax = (*(del[1])) / 2.0;
  293.         }
  294.         else
  295.         if (j == nreg[1] - 1)
  296.         {
  297.             jval = j - 1;
  298.             dy = 1;
  299.             recy = 0;
  300.             ymax = dimv[1];
  301.         }
  302.         else
  303.         {
  304.             start[1] += *(del[1] + j);
  305.             jval = j;
  306.             dy = 0;
  307.             recy = 1.0 / (float) (*(del[1] + j));
  308.             ymax = start[1] + (*(del[1] + j + 1)) / 2.0;
  309.         }
  310.  
  311.         for (y = ymin; y < ymax; y++)    /* y pixels in rgion */
  312.         {
  313. #ifdef REPORT
  314.             if ((loop % 4) == 0)
  315.             {
  316.             printf ("status: %d/192", status);
  317.             fflush (stdout);
  318.             }
  319.             if ((loop++ % 4) == 0)
  320.             sprintf (argv[0], "status:%d/192", status++);
  321. #endif
  322.             cdy = 1.0 - dy;
  323.             start[2] = 0;
  324.             xmin = 0; 
  325.  
  326.             for (k = -1; k < nreg[2]; k++)    /* x regions */
  327.             {
  328.             if (k == -1)
  329.             {
  330.                 kval = 0;
  331.                 dx = 0;
  332.                 recx = 0;
  333.                 xmax = (*(del[2])) / 2.0;
  334.             }
  335.             else
  336.             if (k == nreg[2] - 1)
  337.             {
  338.                 dx = 1;
  339.                 kval = k - 1;
  340.                 recx = 0;
  341.                 xmax = dimv[2];
  342.             }
  343.             else
  344.             {
  345.                 start[2] += *(del[2] + k);
  346.                 kval = k;
  347.                 dx = 0;
  348.                 recx = 1.0 / (float) (*(del[2] + k));
  349.                 xmax = start[2] + (*(del[2] + k + 1)) / 2.0;
  350.             }
  351.  
  352.             /* maps for the interpolation */
  353.             m1 = mapval (jval, kval);
  354.             m2 = mapval (jval, kval + 1);
  355.             m3 = mapval (jval + 1, kval);
  356.             m4 = mapval (jval + 1, kval + 1);
  357.  
  358.             i1 = iptr1 (y, xmin);
  359.             i2 = iptr2 (y, xmin);
  360. #ifdef FOO
  361.             fprintf (stderr,"x, y= %d %d ",xmin,y) ;
  362. #endif
  363.  
  364.             for (x = xmin; x < xmax; x++)    /* x pixels in region */
  365.             {
  366.                 lefty = dy * (m3[*i1]) + cdy * (m1[*i1]);
  367.  
  368.                 righty = dy * (m4[*i1]) + cdy * (m2[*i1]);
  369.  
  370.                 /* CHANGED */
  371.                 tmpint = (int) (righty * dx
  372.                      + lefty * (1.0 - dx));
  373.  
  374.                             if (tmpint > MAXSHORT)
  375.                 {
  376. #ifdef DEBUG
  377.                 fprintf (stderr,"%lf %lf %lf %lf %lf %d\n",lefty,righty,
  378.                      dx,dy,cdy,tmpint) ;
  379. #endif
  380.                                 tmpint = MAXSHORT ;
  381.                 }
  382.                             *i2 = (unsigned char) (tmpint >> 2) ;
  383.                 i1++;
  384.                 i2++;
  385.                 dx += recx;
  386.             }    /* end for x */
  387.             xmin = xmax;
  388.             }        /* end for k */
  389.             dy += recy;
  390.         }        /* end for y */
  391.         ymin = ymax;
  392.         }            /* end for j */
  393. }                /* end */
  394.  
  395. #ifdef godot
  396. report ()
  397. {
  398.     signal (SIGUSR1, SIG_IGN);
  399.     fprintf (stdout, "%3d 192\n", status);
  400.     fflush (stdout);
  401.     signal (SIGUSR1, report);
  402. }
  403. #endif
  404.